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ABSTRACT 

The effect of a newly born star cluster inside a giant molecular cloud (GMC) is 
to produce a hot bubble and a thin, dense shell of interstellar gas and dust swept 
up by the H II expansion, strong stellar winds, and repeated supernova explosions. 
Lying at the inner side of the shell is the photodissociation region (PDR), the 
origin of much of the far- infrared/sub- millimeter /millimeter (FIR/sub-mm/mm) 
radiation from the interstellar medium (ISM). We present a model for the ex- 
panding shell at different stages of its expansion which predict mm/sub- mm and 
far-IR emission line intensities from a series of key molecular and atomic con- 
stituents in the shell. The kinematic properties of the swept-up shell predicted by 
our model are in very good agreement with the measurements of the supershell 
detected in the nearby starburst galaxy M 82. We compare the modeling results 
with the ratio-ratio plots of the FIR/sub-mm/mm line emission in the central 1.0 
kpc region to investigate the mechanism of star forming activity in M 82. Our 
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model has yielded appropriate gas densities, temperatures, and structure scales 
compared to those measured in M 82, and the total H 2 content is compatible with 
the observations. This implies that the neutral ISM of the central star-forming 
region is a product of fragments of the evolving shells. 

Subject headings: galaxies: starburst - galaxies: ISM - ISM: evolution - ISM: 
clouds - ISM: molecules - stars: formation 



1. Introduction 

Starburst is a phenomenon when the star formation rate (SFR) cannot be sustained for 
the lifetime of the galaxy. It is now clear that active star formation or starburst activity is 
common throughout the universe (Madau et al. 1998). The bursts of massive star formation 
can dramatically alter the structure and evolution of their host galaxies by injecting large 
amounts of energy and mass into the ISM via strong stellar winds and repeated supernova 
explosions. The evolution of the superbubbles and supershells that have sizes ranging from 
several tens to hundreds of parsec plays an important role in understanding the amount 
and distribution of warm gas in the ISM. Furthermore, understanding the characteristics of 
starbursts and their relationship with the ISM, as well as to be able to parametrize the star 
formation history are crucial in understanding the galaxy evolution. 

In the past, several models have been used to interpret the infrared, sub-millimeter, 
millimeter line observations of neutral gas in the central regions of nearby starburst galaxies 
(e.g. Mao et al. 2000; Seaquist & Frayer 2000; Wild et al. 1992, and references therein). These 
include the large velocity gradient (LVG) model (Goldreich & Kwan 1974), the steady-state 
PDR model (Tielens & Hollenbach 1985), and the inhomogeneous radiative transfer model 
taking into account non local thermodynamic equilibrium (non-LTE) (Wild et al. 1992). 
These have revealed that the physical conditions (such as gas density, FUV flux, and gas 
kinetic temperature) are enhanced in starburst regions. However, none of these models are 
able to relate the observed line emission properties of the neutral gas in a starburst galaxy to 
its age and star formation history. In this paper, we introduce an evolving starburst model 
for FIR/sub-mm/mm line emission in gas media that allows us to ultimately achieve this 
goal. 

Our model consists of a standard dynamical model of the bubble/shell structure around 
a young star cluster (see Fig. 1), which has been described in many publications (e.g. Castor, 
McCray, & Weaver 1975; Weaver er al. 1977; McCray & Kafatos 1987; Franco et al. 1990; 
Koo & McKee 1992), a time-dependent stellar population synthesis model (Leitherer et al. 
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1999), a fully time-dependent PDR chemistry model (Bell et al. 2005), and a one-dimensional 
non-LTE line radiative transfer model (Rawlings & Yates 2001). In this paper, we conduct 
a preliminary study using this set of models. We first describe the methodology of our 
model (Section 2). We then follow the evolution of a GMC and a swept-up shell induced 
by massive star formation at the center, and calculate the dynamics, thermal structure, and 
the line radiative transfer of the selected molecular and atomic species in the expanding 
shell (Section 3). We compare our modeling results with the observations of the expanding 
supershell and average gas properties in the central 1.0 kpc region of the nearby starburst 
galaxy M 82 (Section 4). Finally, we present the conclusions of this study (Section 5). 

The basic assumptions for our evolving starburst model are (1) star formation occurs 
primarily within the dense optically thick spherical cloud (e.g. Gao et al. 2001), and that all 
stars form instantaneously in a compact spherical cluster located at the center of the cloud 
(the star cluster is therefore treated as a point source), and (2) the starlight produced by the 
central cluster is completely absorbed and reprocessed by the dust in the expanding shell 
(Efstathiou et al. 2000). A summary of our evolving starburst model is presented in Table 1. 



The evolution of a giant molecular cloud is determined by H II expansion in the very 
early stage (t < 10 5 yr), when a hot bubble surrounded by a thin dense shell structure is 
created. The later evolution is driven by the strong stellar winds and repeated supernova 
explosions. We assume that repeated supernova explosions behave like a steady isotropic 
stellar wind injected to the bubble. The hot bubble will eventually cool, and the swept-up 
shell will stall after a few times 10 7 yr. The stars in the young cluster located at the center 
of the GMC are assumed to have masses between 0.1 M and 120 M . The Salpeter initial 
mass function dN/dm* oc m^ 2 - 35 (IMF; Salpeter 1955) is adopted in this study. A top-heavy 
IMF, which has an excess of stars in the mass range 10 - 20 M Q over stars of 5 M or less 
for starburst galaxies (e.g. Rieke et al. 1980), will be investigated in future work. 



The radius and velocity of the H II Expansion due to ionization can be written as 
(Spitzer 1978; Franco et al. 1990), 
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Starburst Models For Gas Media 



2.1. Shell Dynamics 
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where Rs is the initial Stromgren radius in pc, and q ~ 11.5 km s 1 is the sound speed 
in the ionized gas with an equilibrium temperature of ~ 10 4 K. 

Almost as soon as the initial Stromgren sphere is formed, the strong winds start to 
impart large amounts of mechanical energy into the bubble. About 96% of the total wind 
energy is generated by stars with masses > 30 M (McCray & Kafatos 1987). The size of the 
hot bubble is assumed to be much larger than the thickness of the swept-up shell, therefore 
the radius and velocity of the shell in the Winds phase can be written as (McCray & Kafatos 
1987), 

R w (t) = 269.0(^V(t 7 )i, (2a) 
V w (t) = i6.l(^)*(f 7 )-§ (2b) 



where L 38 = L w / (10 38 ergs s^ 1 ), L w is the wind mechanical luminosity, L w = f™ 2 C w C m m2~ 2 ' 35 dm*, 
t 7 = tj (10 7 yr), n is the ambient gas density in cm~ 3 , m 1 = 0.1 M , m 2 = 120 M , C w = 1.0 
x 10 29 , C m = 429.0, and 7 = 3.7 (derived from Abbott 1982). The main-sequence lifetime 
of the most massive star (120 M ) in the star cluster is about 7.0 x 10 5 yr (Mac Low & 
McCray 1988). After this time, we assume that the winds-equivalent energy produced by 
the first supernova and the subsequent ones drives the further expansion of the swept-up 
shell. The radius and velocity of the shell in the Supernova phase can be written as, 



Rs N (t) = 97.0(^) f [(t 7 )f-(^) S ]+^(W), (3a) 
V SN (t) = 5.7(^)W f (3b) 



where iV* is the number of stars with masses > 8 M in the cluster, £51 = E$n / (10 51 ergs 
s _1 ), Esn is the energy produced by each supernova explosion, ti st sN is the time when the 
first supernova occurs in the star cluster, and R w (t lstS N) is the shell radius at ti st sN calculated 
from Equation (2a). The average rate of supernova explosions is ~ 6.3 x 10 35 N*E 5 i ergs s _1 
(McCray & Kafatos 1987). When the energy produced by the stellar winds and/or supernova 
explosions is much greater than the radiative losses, the bubble is adiabatic. This Adiabatic 
phase persists until the radiative cooling becomes important for the hot bubble at t c , 
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t c = 4 x 10 6 Z- L5 (A^£ 51 )iW-^ 



(4) 



where Z is the metallicity with respect to the solar. After t c , the expansion of the 
bubble is no longer energy-driven, but momentum-driven. This momentum-driven phase 
is characterized as the snow-plow (SP) phase. For simplicity, we ignore the momentum 
deposition in the shell by SN ejecta (McCray & Kafatos 1987). Hence, the radius and 
velocity of the shell in the Snow-plow phase can be written as, 



where R c is the radius of the bubble at cooling time t c . The snow-plow phase ends when 
the shell expansion velocity is close to the thermal sound speed of gas in the ISM (typically 
~ 10 km s _1 ). The shell will stall and disperse due to the Rayleigh- Taylor instability (Mac 
Low 1999). 

Our one-dimensional shell dynamical model may overestimate the winds and supernova 
mechanical luminosities as argued recently by Dopita et al. (2005), because the mixing and 
dynamical instabilities will occur in two dimensions, and the ISM is intrinsically inhomo- 
geneous. Dopita et al. (2005) also suggested that the higher ISM pressure in starburst 
regions causes the expanding shell to stall at a smaller radius. Another argument is that 
the gravitational instability may induce new star formation inside the shells. These concerns 
may indicate that the conventional bubble/shell dynamics (Weaver er al. 1977; McCray & 
Kafatos 1987) may need to be modified. 



The PDRs that lie at the inner sides of the clouds or shells centrally illuminated by 
massive star formation are the origin of much of the FIR/sub-mm/mm radiation from the 
ISM. Physical conditions of the swept-up gas in these PDRs are very different from those of 
the cold gas components in the ISM. The gas temperature and density of the swept-up shells 
are a few orders of magnitude higher due to the strong FUV radiation and shock compression. 
The FUV radiation (6 eV < his < 13.6 eV) produced by newborn stars plays an important 
role in the heating and chemistry of PDRs, especially during the early evolution. Other 
sources that may contribute to the shell heating are the mechanical energy input by winds 
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2.2. Physical Conditions of The Swept-up Gas 
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and SN explosions (McKee 1999), shocks caused by the accretion of gas at the outer surface 
of the shell (McKee & Hollenbach 1980), and cosmic rays (Suchkov et al. 1993; Bradford et 
al. 2003). Cosmic rays may play an important role in the heating of swept-up gas after the 
stars with masses > 8 M have terminated as supernovae. Heating sources due to cloud- 
cloud collisions (McCray & Snow 1979) or shell-shell interaction (Scalo & Chappell 1999) 
are not considered in this study. 

The total FUV flux is calculated by integrating the flux of the stellar population spec- 
trum between 912 A and 2055 A for each time step using Starburst99, a time-dependent 
stellar population synthesis model developed by Leitherer et al. (1999). We consider an 
instantaneous burst for the star formation law, where the star formation occurs all at the 
same time (i.e. at age zero). The FUV field strength G incident on the inner surface of the 
shell (visual extinction A v = 0) is then calculated by taking the ratio of the total FUV flux 
to the surface area Anr 2 s (t) of the expanding shell at each time-step. We use the same input 
parameters and assumptions for Starburst99 as those used in the shell dynamics calculation 
(see Table 1). 

The swept-up shell itself is supported by thermal gas pressure and non-thermal pressure 
due to micro-turbulence. The gas temperature decreases toward the outer surface of the 
shell, and the total gas density is assumed uniform. Therefore the pressure is lower at the 
outer surface. The shell density n s refers to the total H_2 density n(H2) in this study. The 
shell density at each time-step is derived from balancing the pressure at the outer surface of 
the shell with the ram pressure, 



„ , A n a v 2 s (t) 

U ° {t) = kT s{ t),, + 5vl (6) 

where n a is the ambient molecular gas density, v s (t) is the expansion velocity, T s (t) is the 
gas temperature at the outer surface of the shell, \x is the mean molecular weight, fi = 0.62 
™>h, m H is the mass of the hydrogen atom, and Svd is the micro-turbulent velocity inside 
the shell. The calculation of the gas temperature profile across the shell will be described 
in the following section. The thickness of the shell d s at each time-step is in turn calculated 
using the continuity equation (or mass conservation law), 
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2.3. The Time-dependent PDR Model 

The gas temperature and chemical abundances of the swept-up shell are calculated self- 
consistently at each depth- and time-step using the time-dependent PDR model developed 
at UCL (called UCL_PDR). A fully time-dependent treatment of the chemistry is employed 
in UCL_PDR which includes 128 species involved in a network of over 1700 reactions (Bell 
et al. 2005, and references therein). The polycylic aromatic hydrocarbons (PAHs) chemistry 
is not included. The reaction rates are taken from the UMIST chemical database (Le Teuff 
et al. 2000). Detailed chemical modeling, heating and cooling mechanisms and the thermal 
balance between them are described in the literature (e.g. Taylor et al. 1993; Papadopoulos 
et al. 2002, and references therein). Heating due to shocks is not included. The UCL_PDR 
code has been modified for the purpose of this study to include a pressure balance check at 
the outer surface of the shell, as well as the evolution information of shell density, thickness, 
and FUV radiation strength. 

The UCL_PDR code assumes a plane-parallel geometry and models the PDR as a semi- 
infinite slab of homogeneous density at a given time-step. The pressure is thus not in 
equilibrium across the PDR region. The FUV radiation field illuminates the shell from one 
side, and it becomes attenuated with increasing visual extinction A v into the shell at a given 
time-step as G = G e~ 2AkAv , where G is the FUV strength at A v = calculated by the 
Starburst99 model. The coefficient 2Ak in front of the A v in the exponent takes into account 
the difference in opacity from the visible to the UV and the influence of grain scattering. 
The timescale for gas in PDRs to reach chemical equilibrium depends on the gas density, 
temperature, degree of ionization, and species involved (Hollenbach & Tielens 1997; van 
Dishoeck & Blake 1998). In our study, this timescale varies from 10 5 to 10 7 yr for the swept- 
up gas. Our comparative tests using a single time-step model fail to reproduce important 
chemical structure features predicted by the fully time-dependent model for ages up to 10 
Myr. The use of a full time-dependent PDR code in which temperature and density changes 
with time is therefore justified in modeling the shell evolution over these timescales. 



2.4. The Non-LTE Line Radiative Transfer Model 

The line radiative transfer properties are calculated using the Spherical Multi-Mol (SM- 
MOL) code. The SMMOL model was also developed at UCL, implementing an accelerated 
A-iteration (ALI) method to solve multi-level non-LTE radiative transfer problems of gas 
inflow and outflow. The code computes the total radiation field and the level populations 
self-consistently. At each radial point, SMMOL generates the level populations and the line 
source functions. Our model assumes that the gas emission originates in the unresolved, 
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homogeneous, spherical expanding shell and that all gas and dust in the H II region have 
been swept up into the shell, i.e. a dustless H II region. The background radiation field 
is assumed to be the cosmic background continuum of 2.73 K. A detailed description of 
the SMMOL radiative transfer model and its implementation can be found in the appendix 
of Rawlings & Yates (2001). The benchmarking and comparison with other line radiative 
transfer models are presented in van Zadelhoff et al. (2002). 

Several programs were developed to separate and extract the gas temperature and frac- 
tional abundances for molecular and atomic species calculated by the UCL_PDR code. These 
extracted gas temperature and abundances, along with the shell density, thickness, radius, 
and expansion velocity computed by the dynamical code, are re-gridded for a spherical geom- 
etry and used as input parameters for the SMMOL code to compute the total line intensity 
or flux. Einstein A and collisional rate coefficients for the molecular and atomic lines are 
taken from the Leiden Atomic and Molecular Database (Schoier et al. 2005). The lowest 10 
energy levels are calculated for all molecules, 3 for atomic [C I] and [O I], and 8 for atomic 
[C II]. 



3. Simulation of An Expanding Shell 

Observational studies have shown that molecular clouds in the Milky Way have a distinct 
mass spectrum Mq MC , with a = -1.5 ± 0.1 (Sanders et al. 1985; Solomon et al. 1987) for 
cloud masses ranging between 10 2 and 10 7 M . Therefore, about 70% of the molecular mass 
in the Galaxy is contained in the GMCs with masses > 10 6 M . These giant molecular 
clouds are known to be associated with active formation of massive stars. If we assume that 
the cloud mass distribution in a starburst galaxy follows a similar index to the galactic one, 
we would expect much of the luminosity of the starburst to arise from the GMCs with a 
fairly narrow range of masses. We adopt a value for the cloud mass of 10 7 M for the GMCs 
in this study. We assume the average gas consumption rate or star-formation efficiency r\ in 
starburst galaxies per 10 8 yr to be 0.25 (Kennicutt 1998). Therefore, the total stellar mass 
M* for the star cluster in the center of the GMC is 2.5 x 10 6 M , and the number of stars 
A* with masses m* > 8 M is about 2.2 x 10 4 . The radius of the GMC is about 50 pc 
with an average cloud density no = 300 cm~ 3 and a cloud core density n c = 2 x 10 3 cm~ 3 
(Plume et al. 1992; Efstathiou et al. 2000). The ambient density n,i sm is assumed to be 30 
cm' 3 (Comeron & Torra 1994). Here we present an idealized case study with this particular 
set of input parameters. 
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3.1. Kinematics of The Swept-up Gas 

The size of the H II region increases slowly with time. The Stromgren radius is about 
4.9 pc assuming the number of Lyman continuum photons is 5 x 10 52 s _1 . The wind bubble 
catches up with the H II ionization front in a time less than 10 5 yr. The strong stellar winds 
cause the bubble to expand quickly into the cloud and sweep up more gas into the shell. 
The total wind power is estimated as L w ~ 1.4 x 10 40 erg s _1 for the star cluster used in 
the model. When the most massive star in the center cluster (120 M ) terminates as a 
supernova at ~ 0.7 Myr, the thin shell caused by the H II region expansion and the stellar 
winds is still expanding at a speed of ~ 40 km s _1 . At this time, the shell has swept up much 
of the mass of its parent cloud, and is propelled into the ISM with a uniform density. The 
mechanical energy produced by the first supernova and the subsequent ones re-energizes the 
shell. A supernova cut-off mass of 8 M Q is assumed. The total energy generated by supernova 
explosions is ~ 2.0 x 10 55 ergs over 40 Myr. At ~ 7.5 Myr, the hot bubble starts to cool 
and loses its internal pressure, at which time the adiabatic phase ends. We adopt 1.0 for the 
metallicity Z with respect to the solar throughout this study. The effect of lower metallicity, 
which is suspected to be present in starburst galaxies, will be discussed in a future paper. 
The radius and velocity of the shell at the end of the adiabatic phase are about 270 pc and 
24 km s^ 1 , respectively. At ~ 50 Myr, the expansion velocity of the shell decreases to ~ 10 
km s _1 , the shell stalls and becomes thicker and less dense. 

Fig. 2 shows the FUV radiation strength Go incident on the inner surface of the shell 
(A v = 0) as a function of time. The Go value is in units of the Habing field, that is 1.6 x 
10~ 3 ergs cm _2 -s~ 1 throughout this study. The FUV strength decreases from about 10 8 to 
10 5 from the onset of star formation to about 5 Myr when most of the massive O stars (> 
30 M ) have terminated as supernovae. It then decreases twice as fast to a value of 40 at 
100 Myr. 

PDRs are the origin of much of the FIR/sub-mm/mm line emission in a starburst 
galaxy. The surface layer (A v ~ 1) contains atomic H, C, C + and O; the transition from 
atomic to molecular hydrogen occurs at the center layer (A v ~ 1 - 2), whilst C + is converted 
into C and then CO over the region A v ~ 2 - 4. H 2 and CO then extend to higher A v 
region and for A v > 10 atomic O begins to be transformed into molecular 2 . The H 2 
molecule provides effective self-shielding from the FUV radiation field. The CO layer also 
shows a degree of self-shielding, and therefore extends deeper into the shell. Small grains 
play an important role in the photoelectric heating of PDRs. Gas heating is dominated by 
collisional deexcitation of FUV-pumped H 2 and vibrationally excited H 2 at the PDR surface. 
The thermal energy radiated by the dust is important for the gas heating at larger optical 
depth (Hollenbach et al. 1991). The gas heating/chemistry at later evolutionary stages is no 
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longer dominated by stellar radiation but by other sources, such as cosmic-rays and X-rays. 
The PDR cooling is dominated by fine structure line emission, such the [C II] 158/xm and [O 
I] 63 /im transitions, whose critical densities are 3 x 10 3 cm -3 and 5 x 10 5 cm -3 , respectively. 
At greater depths, molecular line emission (CO, OH, H 2 0), ro-vibrational transitions of H 2 , 
and gas-dust collisions contribute to the PDR cooling. 

Table 2 summarizes the input parameters for our fully time-dependent PDR model. 
The initial abundance of H 2 is set to n(H 2 )/n# = 0.5 (Hartquist et al. 2003). At the first 
time-step (t — yr) all depth-steps take as their initial abundances the values produced by 
a single-point dense dark-cloud model. The input parameters for the dark-cloud modeling 
are nn = 4 x 10 5 cm -3 , Tqmc = 10 K, Go = 1, and the gas-phase abundances relative to 
H nuclei x He = 7.5 x 10" 2 , x c = 1.8 x 10" 4 , x = 4.4 x 10~ 4 , and x Mg = 5.1 x 10" 6 . 
For subsequent time-steps, the input abundances are re-set to the output abundances of 
the previous time-step generated by the UCL_PDR code. The gas temperature and chemical 
abundances at each depth- and time-step are calculated by balancing the heating and cooling. 
The cosmic-rays ionization rate is enhanced by a factor of 1.5 at later times (t > 10 Myr) to 
artificially include the soft X-rays heating effect on the gas of the shell. We assume that the 
gas-to-dust mass ratio is 100. Fig. 3 shows the shell density n s (or n(H 2 )) and thickness d s as 
a function of time, as calculated by the shell dynamical code and the UCL_PDR code, under 
the condition that the gas pressure at the outer surface of the shell differs from the ambient 
gas pressure by < 10%. The shell density varies between 10 3 and 10 6 cm -3 , and the thickness 
of the shell changes from 10~ 3 to 10 pc over the 100 Myr. We adopt a fixed micro-turbulent 
velocity 5v D = 1.5 km s -1 for the shell. The evolution of the shell density and thickness is 
constrained by the expansion velocity v s , the shell temperature T s , and the ambient density 
n a (See Equation (6) & (7)). Changes in v s and T s are relatively small during the H II 
expansion (n a = Hq or 300 cm -3 ), as a result we see the first plateau as shown in Fig. 3. 
The jump seen at t ~ 2 x 10 4 yr is caused by the change from the H II expansion to the 
Winds phase. During the early Winds phase and before the shell sweeps up all the material 
of its parent GMC (t < 0.8 Myr), the effect due to the shell deceleration is compensated for 
the effect due to the cooling in the shell. This produces a second plateau. After this time, 
the shell expands into a less dense ambient ISM, i.e. n a = n ism or 30 cm" 3 . Less ambient 
pressure causes a decrease in the shell density or a increase in the shell thickness. Fig. 4 
and Fig. 5 show the profiles of the gas temperature and chemical abundances as a function 
of visual extinction A v for an expanding shell at several characteristic ages. The size of the 
PDR changes with time indicated by different maximum values of A v in both Fig. 4 and 
Fig. 5. At ~ 0.7 Myr, all mass in the GMC has been swept into the shell. 
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3.2. 



Molecular and Atomic Line Emission 



The flux and intensity of FIR/sub- mm/mm line emission is calculated for several molec- 
ular and atomic species (CO, HCN, HCO + , C, C + , and 0). The total flux or intensity of each 
line is the sum of the emission from the entire shell. For the initial 0.7 Myr, the emission 
from the parent GMC is also taken into account in the total line emission calculation. In this 
section, we present predictions of the line ratios for CO, [C I], and [C II] for an expanding 
shell. More simulations will be presented and discussed when we illustrate the model by a 
comparison with the observations of M 82 in § 4. 

Molecular CO is known as a good tracer for the diffuse components and total molecular 
gas content in a galaxy, but it is relatively poor tracer of the dense gas directly involved in 
massive star formation. Fig. 6 shows our modeling results for line ratios of high- J transitions 
to the 1-0 transition of the bright and highly abundant 12 CO molecule as a function of the 
starburst age of the shell. 



where is the line intensity, rj+^i is the line intensity ratio, and J = 1,. . .,7. 

For the adiabatic phase (t < 7.5 Myr), strong winds and supernova explosions compress 
the gas in the fast expanding shell to a high density n(H 2 ) > 10 4 cm -3 (see Fig. 3), and the 
strong FUV radiation G > 10 4 heats up the gas and dust of the shell to a temperature > 
100 K (see Fig. 4). A significant amount of highly excited CO line emission is generated 
from the shell and its parent cloud, and therefore the line ratios of r 2 i through are > 1.0. 
At around 10 Myr, all line ratios (1 < J < 7) have dropped below 1.0, the shell has entered 
the snow-plow phase, the corresponding FUV field Go is < 10 4 , the shell density n(H 2 ) is < 
3.0 x 10 3 cm -3 , and the gas temperature in the shell T gas is between 20 and 230 K. 

The far-infrared fine structure lines are the most important cooling lines of the ISM in 
a galaxy. Fig. 7 shows the modeling results of the line intensity ratio of [C II]158/im to [C 
I]610/im and the line flux ratio of [C II]158//m to CO(1-0). About 75% of the [C II]158//m 
emission comes from PDRs and 25% from the H II region (Colbert et al. 1999). The latter 
is not taken into account in our calculations. The [C II]158/im/CO(l-0) line flux ratio rises 
from about 10 to 10 4 after 1 Myr, and then slowly decreases to ~ 10 3 at 80 Myr. It is 
clear that the cooling of the swept-up gas in the expanding shell is dominated by C + , the 
contribution of the CO cooling is a small fraction to the total gas cooling in a massive star 
forming environment. 



rn = h%l ho 



(8a) 

(8b) 
(8c) 
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4. Application to The Nearby Starburst Galaxy M 82 

In section 4.1, we compare our modeling results with the observations of an expanding 
supershell in the nearby starburst galaxy M 82. In section 4.2, we compare our modeling 
results with the average gas properties in the central 1.0 kpc region of this galaxy. 

M 82 is an irregular starburst galaxy located at a distance of about 3.25 Mpc. This 
galaxy has been observed over a wide range of wavelengths. The starburst activity in M 
82 was likely triggered by tidal interaction with its companion M 81 beginning about 10 s 
yr ago in the nucleus, and is currently propagating into the molecular rings. The infrared 
luminosity of M 82 is about 4 x 10 10 L Q arising mostly from the central ~ 400 pc region, 
which has a stellar bar structure and currently has a high supernova rate of ~ 0.05 - 0.1 yr -1 
(Muxlow et al. 1994). The evolutionary scheme in M 82 remains under debate. The most 
common suggested ages of the M 82 starburst in the central regions are 3-7 Myr predicted 
by Colbert et al. (1999) using one instantaneous burst model in dusty media with a 100 
M cut off, and 10 - 30 Myr predicted by Efstathiou et al. (2000) using two instantaneous 
bursts model in dusty media with a 125 M cut off. Recently, Foster-Schreiber et al. (2003) 
presented a more complete evolutionary scheme of the global starburst activity in M 82, and 
suggested that there are two bursts, one occurred at ~ 5 Myr ago and another one at ~ 10 
Myr ago also using instantaneous bursts model in dusty media with a 100 M Q cut off. 



4.1. The Supershell Surrounding The SNR 41.9+58 

Observations have detected an expanding supershell centered around the bright SNR 
41.9+58 in both molecular line and radio continuum (Weiss et al. 1999; Wills et al. 1999). 
This supershell has a diameter of ~ 130 pc, an expansion velocity of ~ 45 km s™ 1 , and a 
mass of ~ 8 x 10 6 M . Using the set of initial cloud conditions selected for our simulation 
(see § 3), i.e. a cloud mass M GMC = 10 7 M , cloud density n = 300 cur 3 , ambient ISM 
density rii sm = 30 cm -3 , and star formation efficiency i] = 0.25, we derive a swept-up shell 
that has very similar characteristics to the observed one. At the observed radius of ~ 65 pc, 
our model indicates an age of 1 Myr, an expansion velocity of ~ 45 km s _1 , and a swept 
up H 2 mass of ~ 7.6 x 10 6 M . The kinetic energy of the observed supershell is estimated 
to be about 1.6 x 10 53 ergs (Weiss et al. 1999). Our model predicts a kinetic energy of ~ 
1.5 x 10 53 ergs for the expanding shell at the age of 1 Myr. The total mechanical energy 
needed for the creation of this supershell is ~ 1.7 x 10 54 ergs, which is contributed by winds 
and supernovae associated with ~ 1700 O stars (> 40 M ). Therefore, about 10% of the 
total energy is present in the form of kinetic energy of the expanding shell. The remarkably 
good agreement between our model results and the observations implies that this supershell 
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may be created by strong winds and supernova explosions from a star cluster with a total 
mass of 2.5 x 10 6 M which occurred in the center about 1 Myr ago. The comparison is 
summarized in Table 3. 

Furthermore, our model predictions of the CO, [C I], and [C II] line ratios for this 
expanding supershell can be used as a comparison with future observations, and also to 
constrain the physical conditions of the gas in the shell (see Fig. 6 and Fig. 7 presented 
in § 3.2). The model line ratios which are greater than 1.0 at t < 8 Myr imply that the 
molecular CO is optically thin in the expanding supershell. Therefore, it is better to look at 
the high CO transitions (J > 3) in this supershell. 

4.2. The Central Starburst Region 

Besides the known expanding supershell centered around SNR 41.9+58, there are other 
undetected shells with sizes ranging from several tens of parsec to more than 1 kiloparsec, 
and kinetic energies ranging from ~ 10 50 ergs to more than 10 54 ergs. These shells would 
likely be present as partial arcs, or fragments, or cloud-like clumps due to strong winds 
and supernova explosions or due to shell-shell and shell-cloud interactions; only a few are 
visible as full circular arcs. The very good agreement between our model and the supershell 
observations indicates that the set of models we have put forward in this paper can be 
used to interpret other shells in a starburst galaxy like M 82. In this section, we illustrate 
the possibilities by comparing our model calculations with the observed FIR/sub-mm/mm 
properties of molecular and atomic line emission in the central starburst regions. 

First of all, Fig. 8 shows the model ratio-ratio plots of HCN(4-3)/(3-2) versus HCN(3- 
2)/(l-0) and HCO+(4-3)/(3-2) versus HCO+(3-2)/(l-0), and a comparison with the observa- 
tions of dense gas in the central 300 pc region (Seaquist & Frayer 2000). The dense gas tracers 
HCN and HCO + are better indicators of active star formation than CO, but poor tracers 
of the total molecular gas content. For both plots, the best agreement is at a starburst age 
of ~ 3 Myr, implying that the expanding shell size is about 300 pc. Secondly, Fig. 9 shows 
the model ratio-ratio plots for 12 CO and 13 CO, and a comparison with the observations of 
the CO line emission from three lobes (north-east, center, and south-west) in the central 300 
to 600 pc regions (Mao et al. 2000). The isotope abundance ratio [ 12 CO]/[ 13 CO] is adopted 
to be 75 for the simulation. The best agreement shown in plots (a), (b), and (d) is at a 
starburst age of ~ 6 Myr, and ~ 3 Myr for plot (c), corresponding to shell sizes between 300 
and 560 pc. The physical conditions for the gas of the shell at age 3-6 Myr are Go ~ 2 x 
10 4 - 15 x 10 4 , n(H 2 ) ~ 1.0 x 10 4 - 2.0 x 10 4 cm" 3 , T gas ~ 50 - 250 K, and total molecular 
gas mass M, mol ~ 0.3 x 10 8 - 2.1 x 10 8 M . Finally, Fig. 10 shows the model ratio-ratio 
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diagram of [0 I]63/xm/[C II]158/im versus [0 I]63/im/[0 I]145/xm, and a comparison with 
the observations of these atomic lines from the central 1.1 kpc region (Negishi et al. 2001). 
The model [C II]158/xm line is underestimated by a factor of about 1.3, since about 25% 
of the total line emission coming from the H II region is not included in the calculation. 
Therefore, the best agreement between the model and the observation is achieved at an age 
of ~ 25 Myr old. The atomic line data are based on a 80" beam whereas the molecular line 
data pertain only to the the lobes and nuclear sources at a 22" beam. Thus, part of the 
reason for the discordant age in the atomic data may be the different regions sampled, since 
they may have a different starburst history. Our predicted gas conditions for the shell at 
this age are G ~ 500, n(H 2 ) ~ 1.8 x 10 3 cm" 3 , T gas > 15 K, and M mol ~ 6.0 x 10 8 M . 
These conditions are consistent with the PDR model fits to the observations by Colbert et al. 
(1999). But the age inferred by Colbert et al. (1999) is 3 - 7 Myr. The large age discrepancy 
between the two different modeling results from the fact that our model includes a more 
massive cluster which then yields the same gas conditions (FUV flux and gas density) at a 
larger distance and hence, in the context of an expanding shell, an older age. It is clear that 
the starburst age of the whole central region is model dependent. More simulations with 
a variety of input cloud conditions and a comparison with data taken at a wider range of 
wavelengths are needed in order to identify the ages of starbursts accurately. 

Although different stages of development are applicable to different central regions of 
M 82, the shell sizes and the physical conditions of the gas within the rings (diameter ~ 
300 - 600 pc) predicted by our model are similar to what is expected from models involving 
expanding shells from a central starburst such as those proposed by Carlstrom & Kronberg 
(1991). Therefore, it is possible that the molecular rings in M 82 are a product of gas that 
was swept-up by the nuclear starburst activity which has evolved for about 10 s yr. Their 
hypothesis is supported by the observations of CO line emission and continuum emission, as 
well as the discovery of supershells that have not yet had time to break out of the galactic 
plane. However, it is important to realize the foregoing interpretation of the lobes as a ring 
or torus is not unique. A number of authors have argued that the molecular rings may be a 
product of Linblad resonance instabilities associated with the gravitational effects of the bar 
(e.g. Shen & Lo 1996; Wills et al. 2000). In future work, we will carry out more simulations 
to test the hypothesis suggested by Carlstrom & Kronberg (1991). 

5. Conclusions 

We have presented a set of starburst models that can be used to relate the observed 
FIR/sub-mm/mm properties of molecular and atomic gas in a starburst galaxy to its age 
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and star formation history. As a preliminary approach, we have illustrated our model by a 
comparison with the observations of the expanding supershell centered around SNR 41.9+58. 
The very good agreement implies that the expanding supershell is created by strong stellar 
winds and SN explosions from a young star cluster (~ 2.5 x 10 6 M ) located in its center. 

Our model predictions of CO, HCN, and HCO + line ratios agree with the molecular 
data for the central lobes (300 - 600 pc) for a shell with an age in the 3-7 Myr range. This 
implies that the molecular rings are possibly a consequence of swept-up or compressed gas 
caused by massive star formation originating in the nucleus of M 82. More simulations in 
future work may be able to justify this hypothesis. The atomic line ratios calculated by our 
model do not fit the observed data as well as the molecular data, but suggest a much older 
shell, because the atomic line emission comes from a much larger region (> 1 kpc). A variety 
of modeling parameters need to be considered to yield more accurate starburst ages. 

Our model also yields appropriate values for the gas density, temperature, and structure 
scales compared to those measured in the center of M 82 (e.g. Lynds & Sandage 1963; Rieu 
et al. 1989; Stutzki et al. 1997; Seaquist & Frayer 2000; Mao et al. 2000; Negishi et al. 2001; 
Ward et al. 2003), and the total H 2 content within the inner 600 pc (~ 2.0 x 10 8 M ) is 
compatible with the observations (e.g. Wild et al. 1992). Therefore, the neutral ISM in the 
central star-forming region of M 82 may be viewed as a product of evolving shells, and is now 
presenting itself in the form of fragments, small cloud clumps, sheets, or even full circular 
arcs. 
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GMC 




Fig. 1. — A schematic structure of an evolving GMC centrally illuminated by a compact 
young star cluster (SC). R s is the radius of the shell, and R b is the radius of the bubble. 
The PDR lies between the thin, dense swept-up shell and the interior (Hollenbach & Tielens 
1997). 
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M GMC : 10 7 M Q , n : 300 cm" 3 , n ism : 30 cm" 3 
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Fig. 2. — Plot of the FUV radiation field strength Go incident on the inner surface of the 
shell (A v = 0) as a function of time (see text for details). 



-21 - 




10 3 10 4 10 5 10 6 10 7 10 8 

Time (yr) 

Fig. 3. — Plot of the shell density (n s or n(H 2 ), solid line) and thickness (d S) dashed line) 
as a function of time for Mqmc = 10 7 M , an initial cloud density ugmc = 300 cm -3 , and 
an ambient ISM density n ism = 30 cm -3 . The radiative cooling of the hot interior occurs at 
t c ~ 7.5 Myr (dotted line). Data for n s and d s after 10 4 years shown in the plots have been 
smoothed. 
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Fig. 4. — The time-dependent gas and dust temperatures as a function of visual extinction 
A v for an expanding shell. Solid lines represent gas temperature, and dashed lines indicate 
dust temperature. 
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Fig. 5. — The time- dependent chemical abundances of the main species (H, H 2 , H + , e~, C, 
C + , 0, and CO) relative to the total hydrogen density, as a function of visual extinction A v 
for an expanding shell. 
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Fig. 6. — Plot of the model CO line intensity ratios of high- J transitions to the 1-0 
transition as a function of starburst age for an expanding shell. The 12 CO line intensities 
are compared in units of K km s -1 . 
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Fig. 7. — Plots of (a) the model line intensity ratio of [C II]158/xm to [C I]610/im as a 
function of time, the line intensities are compared in units of K km s -1 , and (b) the model 
line flux ratio of [C II]158/xm to CO (1-0) as a function of time, the fluxes are compared in 
units of ergs cm~ 2 s _1 . 
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Fig. 8. — The ratio-ratio diagrams of the HCN and HC0 + line intensities. Plots of (a) the 
model HCN(4-3)/(3-2) line ratio versus the HCN(3-2)/(l-0) ratio for a sequence of starburst 
ages: 0, 0.03, 2, 6, 10, 20, 30, 40, 50, and 70 Myr (labeled as 1, 2,. . .,10), (b) the model 
HCO + (4-3)/(3-2) line ratio versus the HCO + (3-2)/(l-0) ratio for a sequence of starburst 
ages: 0.03, 0.7, 2, 4, 10, 20, 30, 40, 50, and 70 Myr (labeled as 1, 2. . .,10). The modeling 
results are indicated by the open circles connected with dotted lines. The filled circles with 
errorbars in the plots are the observed data (see text for details). 
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Fig. 9. — The ratio-ratio diagrams of the 12 CO and 13 CO line intensities. The modeling 
results for a sequence of starburst ages: 0.03, 4, 6, 10, 20, 30, 40, 50, and 70 Myr (labeled as 
1, 2,. . .,9), are indicated by the open circles connected with dotted lines. The filled circles 
connected by solid lines show the observed data for the three lobes in the center of M 82 
(see text for details). 
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Fig. 10. — The ratio-ratio diagram of the fine structure line fluxes. The model [O I]63/xm/[C 
II]158/xm ratio versus [O I]63/im/145/im ratio for a sequence of starburst ages: 0, 0.03, 0.07, 
0.1, 0.7, 4, 6, 8, 10, 20, 40, and 80 Myr (labeled as 1, 2,. . .,12), the line fluxes are compared 
in units of W m -2 . The modeling results are indicated by the open circles connected with 
dotted lines. The filled circles show the observed data for M 82. 



Table 1. Summary of Starburst Models For FIR/sub-mm/mm Line Emission. 



Name 



Description 



Assumptions: 



spherical geometry, non- magnetized (GMCs and shells) 
no interactions between shells, or shell and cloud 
dustless H II regions 

uniform densities of GMCs and ambient media 
all stars form instantaneously, no stars form inside the shells 
stellar mass 0.1 - 120 M with Salpeter IMF dN/dM* oc M~ 2 - 35 
PDRs exist primarily within the expanding shells 



Input Parameters: 



Output Parameters: 



Observational 



GMC mass M GM c = 10 7 M 

average cloud density n = 300 cm~ 3 , cloud core density n c = 2000 cm -3 
ambient ISM density n,i sm = 30 cm~ 3 
star formation efficiency r] = 0.25 
metallicity Z = 1.0 Z 
gas-to-dust ratio = 100 

radius, velocity, temperature, density, and thickness of the shell 
chemical abundances of different molecules and atoms in the shell 
integrated line intensity/flux, line ratios 

line intensities/fluxes and line ratios for molecules (e.g. 12 CO, 13 CO, HCN, HCO + ), 
and atoms (e.g. [C I], [C II], [O I]) 



Note. — See also Table 2 for more input parameters for the time-dependent PDR model. 
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Table 2. Input Parameters For The Time-dependent PDR Model. 



Parameter 



Symbol 



Value 



Starburst age (yr) 

Incident FUV flux (Habing field) 

Turbulent (microturbulence) velocity (km s -1 ) 

PDR surface density (A v = mag) 

Initial gas-phase abundances relative to H a 

PAH abundance 

Dust visual absorption cross section (cm -2 ) 
H 2 formation rate on dust at A v = (cm 3 s -1 ) 
Cosmic-rays ionization rate (s _1 ) 



t 

Go 

5v D 

n H 

xpah 
a v 

c 



o < t < io 8 

10 < G < 10 s 
1.5 

3 < n H < 10 7 



10 



4.0 x 10~ 7 

3.1 x IO" 10 
3.0 x 10~ 18 
1.3 x 10~ 17 



a The initial gas-phase abundances for all depths at the first time-step (t — 
yr) are produced by a single-point dense dark-cloud model (see text for details). 
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Table 3. Characteristics of The Expanding Supershell in M 82. 



Parameter 


Observation 


Model 


Radius (pc) 


65.0 


65.0 


Age (Myr) 


1.0 


1.0 


Expansion velocity (km s^ 1 ) 


45 


45 


Total H 2 molecular gas mass (x 10 6 M Q ) 


8.0 


7.6 


Kinetic Energy (x 10 53 ergs) 


1.6 


1.5 


Total stellar mass in the center cluster (x 10 6 M ) 




2.5 


Total number of stars (> 40 M ) 




1700 


Total Mechanical Energy (x 10 54 ergs) 




1.7 



